home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / skyline / example / dsky.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  4.2 KB  |  172 lines
  1. c
  2. c     This program is a test driver for the unsymmetric skyline solver.
  3. c     This new version does not check the skyline results vs.
  4. c     results from LINPACK.
  5. c     
  6.       integer          ndim, nndim
  7.       parameter        (  ndim = 2500 )
  8.       parameter        ( nndim = ndim*ndim+1 )
  9.       double precision V(nndim), b(ndim), x(ndim)
  10.       double precision normV, normx
  11.       double precision bsave(ndim)
  12.       double precision resid, residn, eps, epsmch, dmach
  13.       double precision ops, sum, time(6), total
  14.       double precision small, dens
  15.       integer          n, job, band, job_setup
  16.       integer          pd(ndim), nr(ndim)
  17.       integer          NumProc
  18.       logical          full
  19.       character*1      afull
  20.       parameter        ( small  = 1.0d-6 )
  21.  
  22.  
  23.       full = .FALSE.
  24.       job = 0
  25.       epsmch = dmach('e')
  26.       eps = 2.d0 * epsmch
  27.  
  28.  
  29. c ___ Read matrix characteristics ___
  30.  
  31.       print *, 'Enter n'
  32.       read( 5, * ) n
  33.       print *, 'Is A full? ( Y/y for yes ):'
  34.       read( 5, 999 ) afull
  35. 999   format(a)
  36.       if( afull .EQ. 'y'  .OR.  afull .EQ. 'Y' ) then
  37.       full = .TRUE.
  38.       else
  39.       print *, 'Enter bandwidth bw in % with respect to n',
  40.      &             'bw <= 0 or bw > 100 generates random band'
  41.       read(5, *) band
  42.       end if
  43.  
  44.  
  45. c ___ Generate matrix and right-hand side ___
  46.          job_setup = 3
  47.          call dsetup( n, V, pd, nr, normV, full, band, b, job_setup )
  48.  
  49.  
  50. c    ___ Save right-hand side for later checking ___
  51.  
  52.      do i = 1, n
  53.         bsave( i ) = b( i )
  54.      end do
  55.  
  56.  
  57. c    ___ Number of floating point operations. ___
  58.  
  59.      if ( full ) then
  60.          ops = (2.0d0*float(n)**3)/3.0d0 + 2.0d0*float(n)**2
  61.      else
  62.         sum = 0
  63.             do i = 1, n
  64.            sum = sum + nr(i)*nr(i)
  65.         end do
  66.         ops = 2.d0*float(sum) + float(pd(n))/2.d0 +
  67.      &          2.d0*float(pd(n))
  68.      end if
  69.  
  70.  
  71. c    ___ Factor the matrix and record the time ___
  72.  
  73.          t1 = second()
  74.       call dskyfa( V, n, pd, info )
  75.      time(1) = second() - t1
  76.  
  77.  
  78.  
  79. c    ___ Solve using the factorization and time it. ___
  80.  
  81.      t1 = second()
  82.       if (job .EQ. 0) call dskysl( V, n, pd, b, 'N', info )
  83.       if (job .NE. 0) call dskysl( V, n, pd, b, 'T', info )
  84.      time(2) = second() - t1
  85.  
  86.  
  87. c    ___ Total time of factorization and solve. ___
  88.  
  89.      total = time(1) + time(2)
  90.  
  91.  
  92. c    ___ Compute a residual to verify results ___
  93.  
  94.      do i = 1, n
  95.         x( i )    = b( i )
  96.          end do
  97.  
  98.          call dsetup( n, V, pd, nr, normV, full, band, b, job_setup )
  99.  
  100.      do i = 1, n
  101.         b( i )  = -bsave( i )
  102.          end do
  103.  
  104.  
  105.      if ( job .EQ. 0 ) then
  106.         do j = 1, n
  107.            jnr = j - nr(j) - 1
  108.            jpd = pd(j) - nr(j) - 1
  109.            do i = 1, nr(j)
  110.               b(j) = b(j) + V( pd(j-1) + i )*x( jnr + i )
  111.            end do
  112.            do i = 1, nr(j)+1
  113.               b(jnr+i) = b(jnr+i) + V( jpd + i )*x( j )
  114.            end do
  115.         end do
  116.     
  117.      else
  118.         do j = 1, n
  119.            jnr = j - nr(j) - 1
  120.            jpd = pd(j) - nr(j) - 1
  121.            do i = 1, nr(j)
  122.               b( jnr+i ) = b( jnr+i ) + V( pd(j-1) + i )*x( j )
  123.            end do
  124.            do i = 1, nr(j)+1
  125.               b( j ) = b( j ) + V( jpd + i )*x( jnr+i )
  126.            end do
  127.         end do
  128.     
  129.      end if
  130.  
  131.  
  132.      resid    = 0.d0
  133.      normx    = 0.d0
  134.      do i = 1, n
  135.         resid = max( resid, abs( b( i ) ) )
  136.         normx = max( normx, abs( x( i ) ) )
  137.          end do
  138.      residn    = resid / ( n * normV * normx * eps )
  139.  
  140.  
  141.      write(6,900)
  142.       write( 6, 901 ) residn, resid, x(1), x(n)
  143. 900      format( / '     norm. resid      resid     ',
  144.      &             '         x(1)          x(n)')
  145. 901      format( 1p4d16.8 )
  146.  
  147.      call skyenvir( 'P', NumProc )
  148.      if ( full ) then
  149.          write(6,902) n
  150.      &                  , NumProc
  151.      else
  152.          dens = pd(n) / ( float(n)**2 )
  153.          write ( 6, 905 ) n,
  154.      &                        NumProc, 
  155.      &                        pd(n), n, dens
  156.      end if
  157. 902      format(// '   Times are reported for matrices of order ', i5
  158.      &           , '  on', i2, ' processor(s)'
  159.      &         )
  160. 905      format(// '   Times are reported for matrices of order ', i5,
  161.      &             '  on', i2, ' processor(s)',
  162.      &        / '       ( density ', i7,'/',i4,'**2  = ', f5.2,' )' )
  163.          write(6,903)
  164. 903      format(/ 6x,'factor',5x,'solve',6x,'total',5x,'mflops',7x)
  165.  
  166.      time(3) = total
  167.      time(4) = ops / (1.d6 * total)
  168.      write(6,904) (time(i), i=1,4)
  169. 904      format( 4(1pe11.3) )
  170.  
  171.      end
  172.